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Abstract 


A new numerical method has been devised and employed to solve a variety of problems 
related to liquid droplet combustion. The basic transport equations of mass, momentum and 
energy have been formulated in terms of generalized nonorthogonal coordinates, which allows 
for adaptive griding and arbitrary particle shape. In this paper example problems are 
solved for internal droplet heating, droplet ignition and high Reynolds number flow over a 
droplet. 


Introduction 


This paper presents initial results of a research effort whose end goals are the 
calculation of single and multiple liquid droplet combustion flows. The complete problem 
of even single droplet combustion presents a severe challenge for present day numerical 
methods because of the multiple space and time scales which can be introduced into the 
problem. These scales are the result of high Reynolds and Peclet numbers as well as flame 
formation around the droplet. Tn order to resolve all the physical phenomena that are 
contained in the problem, it is necessary to use the grid points in a numerical solution 
method very efficiently. A major new advance in the efficient location of grid points is 
the use of generalized nonorthogonal coordinates and adaptive grids. It will be 
shown, and has been shown, that these methods greatly enhance the ability to calculate 
complex flows. 

Tn the present paper major simplifications have been made in the flow equations to 
isolate physical phenomena and test the efficiency of the numerical methods employed. The 
most limiting simplification will be that of constant overall density, however as will be 
shown, the numerical significance of the results are not damaged by the assumptions made. 

Basic transport equations 

The equations for momentum, energy, and mass transport will now be written in terms of 
generalized coordinates, and the numerical methods employed discussed. The starting point 
is the equations in terms of axially-symmetric cylindrical coordinates, and the equations 
for stream function {♦ ) , vorticity (w ) , temperature (T) and reactant species density 

(p ), which are the following: 
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where the following notetion has been used for physical constantsi K 3 -pre*exponential 
reaction constant: ^-activation energy: R-gas constant: h-therauil conduct lvity» p-over- 
all density: and A hm-chesiical heat release. Also^ the independeht variables t# r and i 
are the time, radial coordinate and axial position, respectively. (A more detailed expla- 
nation of the chemical reaction terms will be given later in the paper.) 


A major difficulty with the above formulation is that the surface of the droplet does 
not lie naturally on a constant value of the independent variables r and x. However, this 
difficulty can be re^pved immediately by transforming to generalised nonorthogonal coordi- 
nates C, n, and thereby making the droplet surface correspond to n*0. In order to simplify 
notation it will first be useful to rewrite equations (1) through (S) in a vector form 


)r 9z }r 3s 

where 




Transforming to the variables t , i and n equation (6) becomes 
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where the following new vectors have been defined 


E - j ♦ fep ♦ hj) 

F “ j (6n^ ♦ ♦ TMj) 

R - j ♦ kj) 
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H - ? 
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and the transformation metrics, or areas and volumes, are given by 
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It can easily be seen that the resulting equations are more complex, however the digital 
computer is extremely well equipped to handle this type of algebraic complexity, and with 
some additional programming a much more valuable research tool is obtained. The equations 
as thev are written in (7) will easily handle arbitrary-shaped particles as well as parti- 
cles whose shape is changing as a function of time. However, the major advantage of the 
above formulation is that it allows the use of adaptive grid procedures to be employed, and 
this feature will be shown to be one of the single most important advances in the efficient 
use of numerical methods to solve complex physical problems. 

An interesting feature of these equations is that they allow for the use of variable 
transport properties u, k and D^, although the overall density must be constant. To the 
authors* knowledge this is the first time that a variable transport property formulation of 
the stream f unct ion-vort icity equations has been given. This formulation will be a com- 
plete description of the mass and energy transport processes inside the droplet, whore 
constant overall density is assumed. 


Numerical methods 


The original plan for the numerica) solution method was to employ a fully implicit 
i terat ive scheme to solve tho set of transport equat ions given previously . In order to 
solve these equations an alternatinq-direction- implicit** (ADI) scheme was employed together 
with a Newton procedure to linearize the equations where necessary. The resulting equa- 
tions are block tridiagonal, and the efficient solver of Hindmarsh^ was used to obtain the 
solution. This procedure proved to be unstable because of the large sensitivity that the 
stream function has to all errors. The coupling of the vort icity and stream function equa- 
t ions, plus general ized coordinates, and block solutions along a line, causes the trun- 
cation error to give very inaccurate values of stream function near the boundary. A major 
part of tho problem is due to the fact that the stream function can change by four or five 
orders of magnitude Irom the body surface to the free stream. The values of the stream 
function near the surface are very small and truncation error in the outer part of the flow 
overwhelms the small stream function values near the surface. The solution to the problem 
is point relaxation of tho stream function equation on a previously calculated vorticity 
distribution, followed by iteration between the stream funct ion and vorticity equations. 
This point relaxation method does not couple all the truncation errors together, and very 
good results were obtained , 

The numerical method that was finally employed had the following features: 

^ • S tream fu n ^ i on-vort icity equations 

a. fTt s t -or^er backward, impl icit time differences for time derivatives, 

b. Second-order central differences in space. 

c. ADI solut ion of the vorticity equation, 

d. Point relaxation of the stream function equation. 

e. Global iteration on both equations. 

1 1 . Energ y a nd species equations 

a. The same space and t ime d i f ferencing as the vorticity equation. 

b. Newton linearization of the nonlinear terms. 

c. Block tridiagonal solution with an ADI marching technique. 

The above numerical method is somewhat ad hoc in that it uses two difference methods within 
itself to determine a solution, however the individual methods reflect the numerical nature 
of the equations being solved. The block tridiagonal solution of the energy and species 
equation is excellent for chemical reactions, while the point iteration method allows the 
stream function to converge accurately and smoothly on the vorticity distribution. 

Criterion for grid placement 

The basic criterion for grid placement that was employed in the present paper will now 
be presented. The computational space, l and n, has been normalized so that their numeri- 
cal values go between zero and one, and the grid points are fixed in time. In the physical 
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space the grid points will be placed and moved in time to achieve the resolution of high 
gradient regions. Along a given arc in the physical space the grid points v.«ill be dis- 
tributed in proportion to the gradient of the dependent variable. If the olstance along a 
given arc in physical space is denoted by S, a mathematical statement of th«. relationship 
between the computational and physical space is 

dc . 1^1 a. 

where S is the distance measured on the n « constant arc, and T the dependent variable of 
the transport equation being solved, in this paper the dependent variable is temperature. 
In order to normalize, allow for "optimisation," and remove singularities, the above equa- 
tion is cast into the following form 


C(x,y,t) 


f/s»ax(i + b I— |)d«) 
o * J s * 


where b is an adjustable constant used for "optimisation" of the grid distribution. 

The above equation has some interesting features which will now be discussed briefly. 
For the case b 0 a uniform distribution of points on the nonorthogonal arc is obtained, 
while for large b, constant values of the variable T, or isotherms, are selected. The 
coordinate location equation is solved in an explicit sense at the old time step, and the 
details can be found in the paper by Dwyer et al.^ Also, it should be mentioned that the 
accuracy in solving the equation does not influence directly the accuracy of the finite 
difference solution. With the use of these generalized coordinates and an adaptive grid 
technique a powerful new method is available for numerical solution. 

In the present paper the adaptive griding procedure has only been employed in the 
calculation of the temperature and species concentrations, and has not beci used for the 
stream function-vorticity part of the method. The stream function and vorticity have been 
calculated from a predetermined grid, however the use of adaptive griding for the fluid 
dynamics is being implemented at the present time. 

Results 


The results which will be presented illustrate mainly the capability of the calculation 
method, and do not give a complete description of droplet combustion. The problems which 
have been solved are the following: 

1. Ignition and flame propagation about hot particles. 

2. Separated flow around liquid droplets. 

A description of each of these results will now be given. 

The first problem solved was to calculate the ignition and flame propagation about a 
spherical particle* The reaction mechanism is very simple and consisted of a premixed fuel 
A reacting and going over to species B. The nondimens ional form of the energy and species 
equations are 
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where Pe - 


Peclet Number, based on maximum velocity U. 


inside the droplet 


N0X ” Dimersionless Pre-Exponential Reaction Coefficient. 

6a ” Activation Temperature, Ea/R* 

and it has been assumed that all transport coefficients are constant, and that the thermal 
diffusivity and species diffusivity are equal. 
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For this ignition model problem the velocity field was assumed to be that given by 
Stokes flow, and the overall density was assumed constant. The value for the Peclet number 
was chosen to be 200 and for an initial temperature of premixed reactant of 

The spherical particle surface will be raised impulsively to i^l,0 and the adiabatic 
flame temperature in terms of dimensionless temperature is 

The values for the pre*-exponential coefficient were chosen to be characteristic of 
hydrocarbon fuels, thus the basic time scales of the processes are similar to a moderate 
Reynolds number particle in a premixed hydrocarbon gas (however, the simulation is highly 
simplified as can be easily seen). The first results for the numerical simulation are 
shewn in Figures (1) through ( 4 ) for ■ 2.2*10^, which corresponds to a rather thin 
flame compared to particle diameter. Figures (1) and (2) show the coordinate system (top) 
and isotherm distribution for a uniform grid simulation. A careful study of the isotherm 
distribution shows oscillacions in the flame, and this is due to the high cell Peclet 
numbers. 

The oscillctions are more severe as the ilame moves away from the body, because of the 
natural increase in grid cell sixe occurring in spherical coordinates. As the flame 
approaches the computational boundary there is only one point to describe the i-^aroe struc- 
ture. This lack of resolution results in temperature and species oscillations in the cal- 
culation, incorrect flame speed and eventually termination of tut calculation due to nega- 
tive values of temperature and species. 

A .»econd more accurate calculation with the same number of grid poin^s, but with the 
adaptive grid strategy employed is shown in Figures (3) and (4). The must dramatic feature 
of the calculation is the bunching of the grid points inside the flame structure. This 
removes all the numerical oscillations, and the overall flame speed agrees well with the 
results of Otey and Dwyer. The results of this calculation thus show dramatically the 
usefulness and capabilities of the adaptive grid procedure for combustion problems. 

A second calculation with a reduced reaction rate * 2.2*10^ ) is shown in Figures 

(5) through (8). In this case the flame is much thicker and the coordinate adaption is 
only very slight. However, the ignition and flame propagation processes are much more 
interesting. Figure (5) exhibits the isotherm distribution at an early time with a maximum 
temperature of ^ - 1.0 occur ing at the particle surface. As time increases an ignition 
process occurs at the rear stagnation point (Figure (6)) and then a steady state reaction 
zone is set up on the leeward side of the particle (Figures (7)-(8)). It is easily seen 
from the isotherm distribution that the surface location where the gas temperature rises 
above the wall temperature quickly stabilizes on the leeward side of the particle surface. 
These results show that the present calculation methods will resolve ignition phenomena in 
particle dynamics. 

The final results to be presented are a demonstration of the ability of the numerical 
methods to calculate the fluid flow around and inside droplets. Shown in Figures (9) 
through (11) are the distribution of stream function and vorticity for a solid particle in 
a flow with a Reynolds number of 100, based on particle diameter. Figure (9) exhibits the 
distribution of stream function outside the particle, while Figure (10) shows the stieam 
function Inside the separation bubble (the top figure is the coordinate distributions). A 
more dramatic representation of high Reynolds number Influences can be seen in Figure (11) 
where the normalized vorticity contours are given. The bunching of the contours on the 
windward side of the particle is clear evidence ot the start of boundary layer formation 
and separation. 

All of the above results agree well with the calculations of Le Clair. ^ Figure (12) 
exhibits the streamline pattern inside of a liquid droplet at an external flow Reynolds 
number of 200. For this calculation the ratios of liquid droplet to gas viscosity and 
density are respectively 

lit - • 1000 

M P 

g g 

which are typical values of interest to combustion problems. Therefore, it seems that the 
methods we are employing are ite encouraging, and give strong promise of giving new 
result's for the complete problem of droplet combustion. 

Conclusions 


A r.ew collection of numerical techniques has been assembled to solve problems of heat, 
mass and momentum transport in droplet combustion systems. The major new features of this 
collection of methods are the following: 
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1« Generalized Monorthogonal Coordinates 

2. Adaptive Griding on Temperature Gradients 
3* Block Tridiagonal Solution of Energy and Species Equations 
4, Point Iteration of the Stream Function on the Vorticity Distribution 
All of the above methods have shovn themselves stable and capable of giving improved 
results for the fluid flow# heat transfer, and mass transfer problems solved. 

The physical problems solved in the paper were moderate Reynolds number flow over both 
solid and liquid droplets, as well as a study of ignition around a hot solid particle in a 
Stokes flow. The fluid flow studies reproduced the results of other investigators, and 
thus verified the accuracy of the methods employed. The study of ignition about a hot par- 
ticle showed clearly that Ignition can be delayed until the leeward side of the particle, 
and a flame can be stabilized in the wake of the particle. These results seem to be new, 
and the future inclusion of variable density will allow for a complete description of 
particle ignition. Also, it should be mentioned again that the high-reaction-rate ignition 
studies would not be possible without adaptive griding, because of its efficient use of 
grid points. 


References 

1. Dwyer, H. A., Kee, R. J., Barr, P. K. , and Sanders, B. R. "Transient Droplet Heating 
at High Peclet Number," A.S.M.E. Winter Annual Meeting, Symposium on Computers and Fluids, 
November 1981. 

2. Dwyer, H. A., Kee, R. J., and Sanders, B. R., "Adaptive Grid Method for Problems in 
Fluid Mechanics and Heat Transfer," AIAA J . , Vol. 18, No. 10, October 1980, p. 1205. 

3. Steger, J. L., "Implicit F i n i t e -bl f f e re nee Simulation of Flow About Arbitrary 
Two-Dimensional Geometries," AIAA Journal , vol. 16, July 1978, pp. 679-686. 

4. Roache, P., Computational Fluid Dynamics, Hermosa Publishers, Albuquerque, New 
Mexico, 1971 

5. Hindmarsh, A. C-, "Solution of Block-Tridiagnonal Systems of Linear Algebraic 
Equations," Report UCID-30150, Lawrence Livermore Laboratory, February 1977. 

6. Otey, G. and Dwyer, H. A., "A Numerical Study of the Interaction of Past Chemistry 
and Diffusion," AIAA J. , Vol. 17, June 1979, pp. 606-613. 

7. Le Clair, B. P., Hamielec, A. E,, Pruppacher, H. R. and Hall, W, D., "A Theoretical 
and Experimental Study of the Internal Circulation in Water Drops Palling at Terminal 
Velocity in Air, Journal of Atmoroheric Sciences, V. 29, pp. 728-740. 


36 3 


k 








Ml 


FM 


PostIgnUfQn isotherm distribution 





